/********

This function returns a random number distributed as a Gaussian
on the positive side in the range [mmin,mmax] (or [mzero,mmax] if mmin<mzero) 
with parameters sigma,mzero

exp(-(m-mzero)^2/(2sigma^2))

Based on Normnaldev in Num Rec.
Check explicitly the Gaussian shape with a fit.

*******/


double rangauss(double mzero, double sigma, double mmin, double mmax,
		gsl_rng * rand) {
  

  double u,v,x,y,q,out;
  
  do{    
    //Added one extra exit condition 
    //to sample only the domain [mmin,mmax]
    
    do {
      u = gsl_rng_uniform(rand);
      v = 1.7156*(gsl_rng_uniform(rand)-0.5);
      x = u - 0.449871;
      y = abs(v) + 0.386595;
      q = x*x + y*(0.19600*y-0.25472*x);
    } while (q > 0.27597 && (q > 0.27846 || (v*v) > -4.*log(u)*(u*u)));  
    
    //evaluate the output (take abs to reject only > mzero)
    out = (mzero + abs(sigma*v/u)); 
  } while (out < mmin || out > mmax);
  
  return out;
  
}
